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I.  Introduction 


Several  times  in  science  and  engineering  we  need  to  evaluate  rather  exotic  functions  **  elliptic 
integrals.  Bessel  functions,  and  probability  density  functions,  for  example.  Here's  where  you  can  get 
your  hands  on  some  LISP  code  that  will  evaluate  some  of  these  functions.  The  code  runs  in  both 
MACLISP  and  LISPM  LISP. 

The  functions  in  the  library  arc  scattered  in  several  different  files.  The  descriptions  below  will 
describe  die  function  name,  how  it  should  be  used,  and  the  file  where  that  function  can  be  obtained.  For 
example,  the  function  FAC  TOR  is  in  PS:<GLR.FUNCT>MOD.LSP.  Before  it  can  be  used,  the  code 
must  first  be  loaded  by  doing  one  of: 

(LOAD  " PS : <GLR . FUNC T>MOO . FASL M )  ;  For  MACLISP 

(LOAD  "OZ : PS : <GLR . FUNCT>MOD.QFASL”)  ;  For  MIT  LISPM 

(LOAD  "OZ :  PS :  <GLR .  FUNCT>M00 . Q8W" )  ;  For  Symbolic*  LM l 

(LOAD  "OZ :PS: <GLR . FUMCT>MOD. BIM")  ;  For  Symbolics  3600 

After  die  file  has  been  loaded,  just  call  the  function  with  the  right  arguments: 

(FACTOR  12346670987664321)  *■>  (3333  37  37  333687  333667) 

If  you  have  to  compile  one  of  die  library  files  or  you  want  to  read  the  source  version  of  a  fibraiy 
file  into  your  LISP,  then  you  will  have  to  load  a  special  macro  called  IMPORT-Fllii  that  is  defined  in 
the  file  PS:<Gl.R.LISP>IMPORTJ.SP.  You  may  want  to  use  the  IMPORT-FILE  macro  yourself.  Its 
format  is: 

( IMPORT-FILE  "OZ : PS : <GLR . FUMCT>MOO .,LSP") 

The  argument  is  the  source  filename  with  a  LiSPM  style  hostname  (the  MACLISP  version  of  the  macro 
knows  to  ignore  the  host).  The  behavior  of  the  macro  depends  on  whether  it  is  bdag  compiled, 
evaluated,  or  loaded.  If  the  IMPORT-FILE  form  is  evaluated  or  loaded  (EVAL-WHEN  (EVAL  LOAD) 
...)  and  the  imported  file  has  not  already  been  loaded,  it  will  load  a  binary  file  if  it  can  find  one  or  the 
source  file  if  it  cannot.  At  compile  time  (HVAL-WHKN  (COMPILE) ...),  the  IMPORT-FILE  form  turns 
into 

(PROGN  'COMPILE  <f1r*t  farm  l«  tourco  file*  <lo«4  cod«>) 

By  convention,  the  first  form  in  the  imported  source  file  should  be  a  DECLARE  of  all  die  external 
functions  and  variables  that  the  imported  file  defines.  Including  the  declarations  allows  some  compile 
time  error  checking  and  eliminates  spurious  warnings  about  undefined  functions.  The  <loadcodc> 
conspires  to  load  the  binary  or  source  file  when  when  the  file  we  arc  currently  compiling  gets  loaded. 
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2.  Number  Theory  and  Combinatorial  Functions 

(FACTORIAL  N) 

'The  number  of  ways  to  permute  N  objects. 


PS:<GLR.FUNCT>COMBIN.LSP 


(CHOOSE  N  M)  PS : <GLR .  FUNCT>COMBIN.LSP 

Tlic  number  of  ways  to  place  M  distinguishable  objects  on  N  distinguishable  plates. 

(BELL-NUMBER  N)  ps:<glr.funct>combin.lsp 

The  number  of  ways  to  place  N  distinguishable  objects  on  N  indistinguishable  plates. 

(CATALAN-NUMBER  N)  ps:<gir.funct>combih.lsp 

The  number  of  ways  to  fully  parenthesize  a  string  of  n  symbols. 

(FIB  N)  PS : <GLR . FUNC  T>F IB . LSP 

The  Nth  (F1XNUM)  Fibonacci  number.  FlB(0)=0,  FIB(1)=1 .  Uses  a  log(N)  algorithm. 

(FACTOR  N)  PS:<GLR.FUNCT>MOD.LSP 

Find  the  prime  factors  of  the  integer  N  (FIXNUM  or  BIGNUM).  Returns  a  sorted  list  of  the 

factors. 

(PRIME-TEST  N  ^optional  (TRIALS  60.))  ps:<glr.fumct>mod.lsp 

Tests  the  primality  of  N  using  a  probabilistic  algorithm  (sec  <Solovay>).  Returns  NIL  PRIME, 
orPROBABLY-PRIME(P{crror}  =  2' TR,A,S).  TRIALS  must  be  FIXNUM. 

(SMALLER-PRIME  N)  ps:<glr.fuhct>moo.lsp 

Finds  the  first  prime  that  is  smaller  than  N.  Repeatedly  calls  PRIME-TEST. 


(JACOB I -SYMBOL  P  Q) 

Computes  the  Jacobi  Symbol  of  P  and  Q. 


PS :<GLR. FUNCT>M0D . LSP 


(TOTIENT  N)  r 

Computes  Euler’s  totient  function  of  N.  Calls  FACTOR  as  a  subroutine. 


PS:<GLR. FUNCT>MOD.LSP 


3.  Discrete  Fourier  Transform 

These  functions  arc  used  to  compute  discrete  fouricr  transforms  <Oppcnhcimcr>.  Ihc  inputs 
arc  two  FLONUM  arrays  that  represent  the  real  and  imaginary  parts  of  the  input  The  algorithms  require 
input  arrays  for  the  transform  to  be  bit  reversed;  there  are  functions  for  doing  the  reversal.  I  he  functions 
also  require  some  initialization.  'Ihc  length  of  the  DFT  must  be  a  power  of  2. 
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3  DISCKirn-  IOURII-R  TRANSl'ORM 


X[k]  =  2  ^  cxp(- j  2  irkm/N) 
x{k)  =  (1/N)  2  N=1q  Cxp(  j  2  it  k  m  /  N) 

(DFT-INIT  LOGN)  PS:<GLR.FU*CT>OfT.l.SP 

DFT-INIT  returns  a  structure  containing  some  tables  that  arc  needed  by  the  DFT  algorithm. 
This  structure  should  be  saved  away  and  given  to  the  DFT  routine  each  time  it  is  called.  The  structure 
need  only  be  computed  once  for  each  size  DPT.  The  length  of  the  DFT  is  2*®^^. 

(DFT-FORWARD  X-REAL  X-IMAG  TABLES)  ps:<glr.fumct>oft.lsp 

T  his  function  docs  a  discrete  Fourier  transform.  The  results  arc  written  back  into  X-REAL  and 
X-IMAG  (the  previous  contents  arc  lost)  and  arc  in  sequential  order  (ic  --  not  bit  reversed).  TABi.ES  is 
as  returned  by  DFT-INIT. 

(DFT-REVERSE  X-REAL  X-IMAG  TABLES)  ps:<glr.funct>oft.lsp 

T  his  is  just  like  DFT-FORWARD.  but  it  docs  the  inverse  transform.  The  input  X  arrays  should 
be  in  sequential  (not  bit-reversed)  order.  The  resulting  arrays  arc  also  in  sequential  order. 

(DFT-REVERSE -ARRAY  ARRAY  TABLES)  ps:<6lr.fwict>oft.lsp 

Ibis  bit  reverses  the  first  21 elements  of  ARRAY. 

(DFT-610  X-REAL  X-IMAG  TABLES)  ps:<glr.fumct>oft.lsp 

Ibis  function  docs  a  forward  transform.  X-REAL  and  X-IMAG  are  bit-reversed  input  arrays 
(to  reverse  them,  use  DFT-RKVERSB-ARRAY).  'Ibe  results  are  returned  in  X-REAL  and  X-IMAG  in 
sequential  order. 

(DFT-618  X-REAL  X-IMAG  TABLES)  ps:<glr.funct>oft.lsp 

Ibe  reverse  transform.  Input  arrays  are  sequential  order,  output  arrays  are  bit-reversed.  Use 
DFT-REVERSE-ARRAY  to  put  them  in  sequential  order. 

4.  Trigonometric  and  other  Functions 

These  approximations  come  from  <Abramowitz>,  <DEC>,  and  <Hastings>. 

(EXP10  X)  PS:  <GLR.  FUNCT>EXTFCN.LSP 

Computes  lOtX. 


4  TRIGONOMETRIC  AND  OTHER  JUNCTIONS 
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(LOGIO  X).  PS : <GLR . FUNCT>EXTFCN . LSP 

Log  base  lOofX. 

(TAN  X) 

(SEC  X) 

(CSC  X) 

(COT  X) 

(ASIN  X) 

(ACOS  X) 

(ASEC  X) 

(ACSC  X) 

(ACOT  Y  X)  PS:<GLR. FUNCT>HYPER .LSP 

Various  trigonometric  functions. 

(ERROR-FUNCTION  X)  PS : <GLR . FUNCT>EXTFCN . LSP 

Error  function.  eps<  1.5E-7. 

(2/sq  rt(  ir ))  /0xcxp(-t2)dt 

(BESSEL-I  N  X)  PS:<GLR.FUNCT>BESSEL.LSP 

Modified  Bessel  function  for  integer  order  N.  eps  <  1.6E-7 

(BESSEL-J  N  X)  PS:<6LR.FURCT>BESSEL.LSP 

Bessel  function  for  integer  order  N.  eps  <  5E-8. 

(GAMMA- FUNCTION  A)  PS : <GLR . FUNCT> GAMMA . LSP 

Gamma  function. 

r<a)=/0°°lr'c>dt 

(GAMMA- FUNCTION- INCOMPLETE  A  X)  PS : <GLR . FURCT>GAMMA . LSP 

Incomplete  Gamma  function.  0  <  =  x  <  =  inf,  but  breaks  if  X  is  taigc  because  of  roundoff  error. 

T(a,x)  =  /QXta*le‘tdt 

(BETA-FUNCTION  A  B)  PS : <GLR . FUMCT>GAMMA . LSP 

Beta  function.  A  and  B  arc  FLONUM. 

a.  b)  =  /q1  ta‘l  (l-t)b  l  dt ,  a>0.  b>0 

(BETA-FUNCTION-INCOMPLETE  A  B  X)  P$:<6LR. FUNCT>GAMMA.LSP 

Incomplete  Beta  function.  A  and  B  arc  FLONUM,  but  one  of  A  or  B  must  be  an  integer  (eg, 
3.0). 

b. x)  =  /0* t*'1  (l-t)b*1  dt,  0<a, 0<b,(K=x<=l 
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4  TRIGONOMI-TRIC  ANDOT1IUR IOJ  NOTIONS 


(ELLIPTIC-INTEGRAL-K  M)  ps^glr.fuuct^llip.lsp 

Flliptic  integral  of  the  first  kind.  k(M).  0<=M<1.  cps  =  2.0E-8. 

/0W/2  (1  -  m  sin2(0))'°-5  d 9 

(ELLIPTIC-INTEGRAL-KC  M)  PS :  <GLR .  FUNCT>ELLIP .  LSP 

Complementary  elliptic  integral  of  the  first  kind,  k’(M).  k’(M)  =  k(l-M). 

( ELLIPTIC- INTE6RAL-E  M)  ps:<glr.fumct>ellip.lsp 

Flliptic  integral  of  the  second  kind,  c(M).  cps  <  2.0E-8. 

/  sqrt(  1  * m  sm2(0))  &0 

(ELLIPTIC-INTE6RAL-EC  M)  PS : <GLR . FUHCT>ELLIP . LSP 

Complementary  elliptic  integral  of  the  second  kind,  e*(M).  e’(M)  =  c(l-M). 

(ELLIPTIC-SINE  u  M)  PS : <GLR . FUNCT>ELLIP . LSP 

elliptic  sine  function  (SN(u,M» 

(ELLIPTIC-COSINE  u  M) )  p$:<glr.furct»ellip.lsp 

elliptic  cosine  function  (CN(u,M)) 


5.  Linear  Regression 

Say  we  arc  trying  to  fit  some  (x,  y)  data  to  a  function  that  looks  like: 

V (x)  -  aO*  1  ♦  al'gl(x)  ♦  *2*g2(x)  ♦  ...  ♦  an*gn(x) 

where  the  a(i)  arc  constants  that  we  are  trying  to  determine,  the  gpl  are  linearly  independent  functions 
that  we  specify,  and  YTx)  is  the  value  of  the  fitted  equation.  Such  a  fit  is  called  a  linear  regression  and 
here  arc  some  functions  that  will  do  it 

(FIT  FCTN  X  Y  N  (optional  (MODE  0))  ps:<glr.fumct>fit.lsp 

Let  X  be  a  one  dimensional  array  of  x  values  (the  x  values  could  themselves  be  vectors).  Y  is  a 
one  dimensional  array  of  FLONUM  values  for  the  corresponding  X  input  N  is  the  number  of  terms  in 
the  equation  we  arc  fitting  (ic,  n  in  the  equation  above).  FCTN  is  a  function  of  N+l  arguments  that 
evaluates  the  equation  above.  That  is  FCTN  looks  like: 

(LAMBDA  (X  AO  At  A2  AS  . . .  A») 

...  ) 

and  calculates  the  above  equation.  FIT  returns  a  Ksi  of  the  coefficients  Ai  as  the  COR  of  its  value. 

FIT  actually  uses  the  procedure  Rl-GRFSS  in  PS:<GI  J<.FUNCI>RHGRF5.I5P,  which  is 
modeled  after  CBevingtonX  If  you  want  to  do  something  a  little  more  complicated,  then  ask  me  about 
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that  procedure. 


6.  Functional  Minimization 

( FHFP  FUNCT  N  X  G  EST  EPS  LIMIT)  PS :<GLR. FUNCT>FMFP. LSP 

Fletchcr-Powell's  functional  minimization  procedure  (adapted  from  <Kucstcr>  and 
<l.ucnbcrgcr>).  This  finds  a  minimum  of  a  multivariate,  unconstrained,  nonlinear  function.  That  is,  find 
the  vector  X  such  that  F(X)  is  a  local  minimum,  'fhc  arguments  to  this  function  arc: 

N  number  of  independent  variables 

X  vectorfN]  of  initial  variable  values 

contains  the  result  vector  on  exit 
G  vector[N]  in  which  to  store  the  gradient 

EST  estimate  of  minimum  value  of  objective  function 
EPS  test  value  representing  the  expected  absolute  error  In  movement 
LIMIT  maximum  number  of  iterations 

FUNCT  user  supplied  objective  function  that  computes  F  and  the  gradient 
(FUNCT  N  X  G)  returns  F(X)  and 
also  fills  in  the  vector[N]  G  with  the  gradient 

(MARQUARDT  NKXYZ  FUNC  DERIV  B  BMIN  BMAX  BV) 

PS : <GLR . FUNCT>MARQ . LSP 

Marquardf s  parameter  fitting  procedure  (adapted  from  <Kuestcr>).  Given  N  data  points  (X|iJ, 
Y[i]>  and  a  function  F  with  K  parameters  U(j],  find  the  K  parameters  B[j]  that  best  fit  die  data,  'fhc  fitted 
values  are  returned  in  array  Z.  The  arguments  are: 

N  —  number  of  data  points  (FIXNUM) 

K  —  number  of  unknowns  (FIXNUM) 

B  —  vector[K]  of  (FLONUM)  unknowns 

BMIN  --  vector[KJ  of  (FLONUM)  minimum  values  of  B[j] 

BMAX  —  vector[K]  of  (FLONUM)  maximum  values  of  B[j] 

X  —  vector[N]  Independent  variable  data  points 
(this  vector  might  be  a  vector  of  vectors) 

Y  —  vector[N]  of  (FLONUM)  dependent  variable 

l  --  vectorfN]  of  (FLONUM)  computed  values  of  dependent  variable 

BV  --  vector[K]  of  (FLONUM)  codes 

1.0  ->  numerical  derivatives 
0.0  ->  do  not  change  this  unknown 
FUNC  —  (FUNC  X  K  B  N  Z  ZOFF)  computes  the  function  F  (using  the 
X  parameters  In  vector  B)  for  each  of  the  X[0]  ...  X[N-1] 
and  puts  the  respective  results  Into  Z[ZOFF+0]  ...  Z[Z0FF+N-1] 

0ERIV--  (not  used) 
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7 IIYPI  RHOUC  &  INVI-RSI:  IIYPURUOI JC  FUNCTIONS 


7.  I  lypcrbolic  &  Inverse  I  lyperbolic  Functions 

lliis  functions  use  the  Common  Lisp  <Stccic>  names,  but  arc  only  defined  for  real  arguments. 
They  call  KXP  to  compute  their  results. 

(COSH  X) 

(SINH  X) 

(TANH  X) 

(COTH  X) 

(SECH  X) 

(CSCH  X)  PS :  <GLR .  FUNCT>HYPER. LSP 

Hyperbolic  functions.  FLONUM  argument  and  FLONUM  result 

(ACOSH  X) 

(ASINH  X) 

(ATANH  X) 

(ACOTH  X) 

(ASECH  X) 

(ACSCH  X)  PS :  <GLR .  FURCT>HYP'R .  tSP 

Inverse  hyperbolic  functions.  FLONUM  argument  and  FLONUM  result 

8.  Numerical  Integration  j 

(INTEGRATE-TRAPEZOIDAL  F  XO  XI  N)  ps : <glr . funct> inte&r . lsp 

Uses  die  trapezoidal  iuIc  to  integrate  the  function  F  (of  one  FLONUM  argument)  from  XO  to 
XI  with  N  (I-IXNUM)  iterations. 

(INTEGRATE-SIMPSON  F  XO  XI  N)  PS : <GLR . FUNCT>INTE6R . LSP 

Uses  Simpson’s  rule  to  integrate  the  function  F  (of  one  FLONUM  argument)  from  XO  to  XI 
with  N  (F1XNUM)  iterations. 

(INTEGRATE-EULER  F  XO  YO  H  XI)  ps:<glr.furct>rurge.lsp 

Uses  the  forward  Killer  method  to  integrate  F(x,  y)  from  position  (XO,  YO)  to  a  position  (XI,  Yl) 
using  a  step  size  of  H.  Yl  is  the  value  of  1NTKGRATE-EUI.HR.  F  is  a  function  ofX  and  Y  (FLONUM) 
and  should  return  DY/DX  for  die  given  position. 

(INTEGRATE-RUNGE-KUTTA  F  XO  YO  H  XI)  ps:<glr.funct>runge.lsp 

This  is  just  like  KUI.F.R  except  it  uses  a  fourth  order  RUNGE-KUTI'A  method  instead  of 
Filler's  method. 


8  NUMERICAL  INTEGRATION 
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( INTEGRATE-MULTI-STEP-H  F  XO  YO  H  XI 

&optiona1  (CD  1.0E-8)  (CH  1.0E-6))  PS : <GLR . FUNCT>RUNGE . LSP 

Uses  a  multi-step  predictor-corrector  method  to  integrate  DKRIV  (sec  <Hamming>  page  407). 
Integrates  from  XO  to  XI  using  an  initial  step  size  of  H.  When  the  routine  estimates  that  the  error  is 
below  CD,  the  step  size  is  doubled:  when  the  error  is  above  CH,  the  step  size  is  halved.  For 
reasonableness,  CD  «  CH. 

9.  Root  Finding  Functions 

(BISECTION-ROOT-SEARCH  F  XO  XI)  ps:<glr.funct>root.l$p 

Uses  bisection  search  to  find  a  zero  of  F(X)  between  XO  and  XI. 

( FALSE-POSITION-SEARCH  F  XO  XI  EPS)  PS : <GLR . FUNCT>ROOT . LSP 

Uses  false  position  search  to  find  a  zero  of  F(X)  with  initial  guesses  XO  and  XL 

(FALSE-POSITION-CONVERGE  F  XO  XI  EPS)  PS : <GLR . FUNCT>R00T . LSP 

Uses  false  position  convergence  to  find  an  X  =  F(X)  with  initial  guesses  XO  and  XI. 

(CONVERGE  F  XO  EPS)  pS:<glr.fdnct>root.lsp 

Find  an  X  =  F(X)  with  initial  guess  of  XO.  Uses  Wcgstcin’s  method. 

10.  Probability  and  Statistics 

The  following  functions  arc  useful  in  probability  and  statistics.  There  arc  functions  for 
computing  probability  density  functions  p(x),  computing  cumulative  distributions  P(x).  and  generating 
random  numbers  from  a  particular  distribution.  Most  of  the  approximations  come  from  <Abramowitz>. 
Most  of  the  generators  use  the  ratio  of  uniform  deviates  method  <Kindcrman>. 

(UNIFORM-DENSITY  X) 

(UNIFORM-CUMULATIVE  X) 

(UNIFORM-RANDOM-NUMBER)  PS: <6LR . FUNCT>STATIS. LSP 

Probability  functions  for  die  uniform  distribution. 

p(x)  =  l  0<=  x  <=  1 

(NORMAL-DENSITY  X) 

(NORMAL-CUMULATIVE  X)  ps:<6lr.furct>stati$.lsp 

Probability  functions  for  normal  (Gaussian)  distribution. 

p(x)  =  (l/sqrt(2  w))  exp  (*x^/2) 
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(NORMAL-TAIL  ALPHA)  ps:<glr.fu*ct>stat!S.lsp 

Given  a  probability  ALPHA,  find  the  X  such  that  ALPHA  =  1-P(X),  where  P(X)  is  the 
cumulative  distribution,  t  his  function  provides  a  one-sided  tail. 

(NORMAL-RANDOM-NUMBER)  PS:<GLR.Fui«CT>STAns.LSP 

Generates  a  random  number  from  the  unit  normal  distribution. 

(EXPONENTIAL-DENSITY  X  LAMBDA) 

( EXPONENTIAL-CUMULATIVE  X  LAMBDA) 

(EXPONENTIAL-RANDOM-NUMBER  LAMBDA)  ps:<glr.funct»statis.lsp 

Exponential  probability  functions. 

p(x.  X)  =  X  exp(-  X  x),  X  >  0 

( POISSON-DENSITY  N  TIME) 

(POISSON-CUMULATIVE  N  TIME) 

( POISSON-RAfJDOM-NUMBER  TIME)  PS : <GLR . FUMCT>STATIS . LSP 

Poisson  probability  functions.  Probability  of  exactly  N  (FfXNUM)  arrivals  within  TIME 
(FLONUM)  for  a  Poisson  process  with  an  average  arrival  rate  (X)  of  1. 

p(n.  ra)  =  (mn  cxp(-m))/nl,  m~  X  t 

(CHI-SQUARE-DENSITY  X  N) 

(CHI-SQUARE-CUMULATIVE  X  N) 

(CHI-SQUARE-RANDOM-NUMBER  N)  PS : <GLR . FUNCT>STATIS . LSP 

The  CHI-SQUARE  distribution  with  N  (FIXNUM)  degrees  of  freedom. 

(T-Dr  :3ITY  X  N) 

(T-CUMULATIVE  X  N) 

(T-RANDOM-NUMBER  N)  PS : <GLR . FUNCT>STATIS . LSP 

Student's  T  distribution  with  N  (FIXNUM)  degrees  of  freedom.  N  must  be  even  in 
T-CUMULATIVE. 

(T-TWO-SIDED  X  N)  PS : <GLR . FUNCT>STATIS . LSP 

For  the  T  distribution,  the  probability  that  T  falls  within  -X  to  +X.  N  must  be  an  even 
FIXNUM. 

(F-DENSITY  X  M  N) 

(F-CUMULATIVE  X  M  N) 

(F -RANDOM-NUMBER  M  N)  PS : <GLR . FUNCT>STATIS . LSP 

Sncdccors  F  distribution  with  M  and  N  (FIXNUM)  degrees  of  freedom. 
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II 


(GAMMA-DENSITY  X  A) 

(GAMMA-CUMULATIVE  X  A) 

(GAMMA-RANDOM-NUMBER  A)  PS : <GLR . FUNCT>STATIS . LSP 

Gamma  probability  functions  with  parameter  A  (Fl.ONUM).  ’Hie  random  number  generator 
must  have  A  >  1.0. 

p(x,  a)  =  (l/r(a))  x3'^  c'x,  a>0 

(BETA-DENSITY  X  A  B) 

(BETA-CUMULATIVE  X  A  B) 

(BETA-RANDOM-NUMBER  A  B)  PS . «GLR . FUNCT>STATIS. LSP 

Beta  probability  functions  with  parameters  A  and  B  (FLONUM).  The  random  number 
generator  is  slow,  so  keep  A  and  B  small  (say  below  10). 

p(x,  a,  b)  =  (l//?(a,b»  x3"1  (l-x)b*1,  a>0.  b>0 

(CAUCHY-DENSITY  X) 

(CAUCHY-CUMULATIVE  X) 

(CAUCHY-RANDOM-NUMBER)  ps.<glr.funct>statis.lsp 

Probability  functions  for  the  Cauchy  distribution. 

p(x)=  l/(7T(l  +  X2)) 
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